home *** CD-ROM | disk | FTP | other *** search
/ Aminet 31 / Aminet 31 (1999)(Schatztruhe)[!][Jun 1999].iso / Aminet / dev / c / random.lha / random / r250 / r250.c < prev    next >
C/C++ Source or Header  |  1999-02-28  |  8KB  |  389 lines
  1. /******************************************************************************
  2.  
  3. *  Module:  r250.c   
  4.  
  5. *  Description: implements R250 random number generator, from S. 
  6.  
  7. *  Kirkpatrick and E.  Stoll, Journal of Computational Physics, 40, p. 
  8.  
  9. *  517 (1981).
  10.  
  11. *  Written by:    W. L. Maier
  12.  
  13. *
  14.  
  15. *    METHOD...
  16.  
  17. *        16 parallel copies of a linear shift register with
  18.  
  19. *        period 2^250 - 1.  FAR longer period than the usual
  20.  
  21. *        linear congruent generator, and commonly faster as
  22.  
  23. *        well.  (For details see the above paper, and the
  24.  
  25. *        article in DDJ referenced below.)
  26.  
  27. *
  28.  
  29. *    HISTORY...
  30.  
  31. *        Sep 92: Number returned by dr250() is in range [0,1) instead 
  32.  
  33. *            of [0,1], so for example a random angle in the
  34.  
  35. *            interval [0, 2*PI) can be calculated
  36.  
  37. *            conveniently.  (J. R. Van Zandt <jrv@mbunix.mitre.org>)
  38.  
  39. *        Aug 92: Initialization is optional.  Default condition is 
  40.  
  41. *            equivalent to initializing with the seed 12345,
  42.  
  43. *            so that the sequence of random numbers begins:
  44.  
  45. *            1173, 53403, 52352, 35341...  (9996 numbers
  46.  
  47. *            skipped) ...57769, 14511, 46930, 11942, 7978,
  48.  
  49. *            56163, 46506, 45768, 21162, 43113...  Using ^=
  50.  
  51. *            operator rather than two separate statements. 
  52.  
  53. *            Initializing with own linear congruent
  54.  
  55. *            pseudorandom number generator for portability. 
  56.  
  57. *            Function prototypes moved to a header file. 
  58.  
  59. *            Implemented r250n() to generate numbers
  60.  
  61. *            uniformly distributed in a specific range
  62.  
  63. *            [0,n), because the common expedient rand()%n is
  64.  
  65. *            incorrect.  (J. R. Van Zandt <jrv@mbunix.mitre.org>)
  66.  
  67. *        May 91: Published by W. L. Maier, "A Fast Pseudo Random Number 
  68.  
  69. *            Generator", Dr. Dobb's Journal #176.
  70.  
  71. ******************************************************************************/
  72.  
  73.  
  74.  
  75. #include <stdlib.h>
  76.  
  77. #include "r250.h"
  78.  
  79.  
  80.  
  81. /**** Static variables ****/
  82.  
  83. static int r250_index = 0;
  84.  
  85. static unsigned int r250_buffer[250] = {
  86.  
  87.     15301,57764,10921,56345,19316,43154,54727,49252,32360,49582,
  88.  
  89.     26124,25833,34404,11030,26232,13965,16051,63635,55860,5184,
  90.  
  91.     15931,39782,16845,11371,38624,10328,9139,1684,48668,59388,
  92.  
  93.     13297,1364,56028,15687,63279,27771,5277,44628,31973,46977,
  94.  
  95.     16327,23408,36065,52272,33610,61549,58364,3472,21367,56357,
  96.  
  97.     56345,54035,7712,55884,39774,10241,50164,47995,1718,46887,
  98.  
  99.     47892,6010,29575,54972,30458,21966,54449,10387,4492,644,
  100.  
  101.     57031,41607,61820,54588,40849,54052,59875,43128,50370,44691,
  102.  
  103.     286,12071,3574,61384,15592,45677,9711,23022,35256,45493,
  104.  
  105.     48913,146,9053,5881,36635,43280,53464,8529,34344,64955,
  106.  
  107.     38266,12730,101,16208,12607,58921,22036,8221,31337,11984,
  108.  
  109.     20290,26734,19552,48,31940,43448,34762,53344,60664,12809,
  110.  
  111.     57318,17436,44730,19375,30,17425,14117,5416,23853,55783,
  112.  
  113.     57995,32074,26526,2192,11447,11,53446,35152,64610,64883,
  114.  
  115.     26899,25357,7667,3577,39414,51161,4,58427,57342,58557,
  116.  
  117.     53233,1066,29237,36808,19370,17493,37568,3,61468,38876,
  118.  
  119.     17586,64937,21716,56472,58160,44955,55221,63880,1,32200,
  120.  
  121.     62066,22911,24090,10438,40783,36364,14999,2489,43284,9898,
  122.  
  123.     39612,9245,593,34857,41054,30162,65497,53340,27209,45417,
  124.  
  125.     37497,4612,58397,52910,56313,62716,22377,40310,15190,34471,
  126.  
  127.     64005,18090,11326,50839,62901,59284,5580,15231,9467,13161,
  128.  
  129.     58500,7259,317,50968,2962,23006,32280,6994,18751,5148,
  130.  
  131.     52739,49370,51892,18552,52264,54031,2804,17360,1919,19639,
  132.  
  133.     2323,9448,43821,11022,45500,31509,49180,35598,38883,19754,
  134.  
  135.     987,11521,55494,38056,20664,2629,50986,31009,54043,59743
  136.  
  137.     };
  138.  
  139.  
  140.  
  141. static unsigned myrand();
  142.  
  143. static void mysrand(unsigned newseed);
  144.  
  145.  
  146.  
  147. /**** Function: r250_init  
  148.  
  149.     Description: initializes r250 random number generator. ****/
  150.  
  151. void r250_init(int seed)
  152.  
  153. {
  154.  
  155. /*---------------------------------------------------------------------------*/
  156.  
  157.     int        j, k;
  158.  
  159.     unsigned int mask;
  160.  
  161.     unsigned int msb;
  162.  
  163. /*---------------------------------------------------------------------------*/
  164.  
  165.     mysrand(seed);
  166.  
  167.     r250_index = 0;
  168.  
  169.     for (j = 0; j < 250; j++)     /* Fill the r250 buffer with 15-bit values */
  170.  
  171.         r250_buffer[j] = myrand();
  172.  
  173.     for (j = 0; j < 250; j++)     /* Set some of the MS bits to 1 */
  174.  
  175.         if (myrand() > 16384)
  176.  
  177.             r250_buffer[j] |= 0x8000;
  178.  
  179.     msb = 0x8000;       /* To turn on the diagonal bit   */
  180.  
  181.     mask = 0xffff;      /* To turn off the leftmost bits */
  182.  
  183.     for (j = 0; j < 16; j++)
  184.  
  185.         {
  186.  
  187.         k = 11 * j + 3;             /* Select a word to operate on        */
  188.  
  189.         r250_buffer[k] &= mask;     /* Turn off bits left of the diagonal */
  190.  
  191.         r250_buffer[k] |= msb;      /* Turn on the diagonal bit           */
  192.  
  193.         mask >>= 1;
  194.  
  195.         msb >>= 1;
  196.  
  197.         }
  198.  
  199. }
  200.  
  201.  
  202.  
  203. /**** Function: r250 Description: returns a random unsigned integer k
  204.  
  205.                 uniformly distributed in the interval 0 <= k < 65536.  ****/
  206.  
  207. unsigned int r250()
  208.  
  209. {
  210.  
  211. /*---------------------------------------------------------------------------*/
  212.  
  213.     register int    j;
  214.  
  215.     register unsigned int new_rand;
  216.  
  217. /*---------------------------------------------------------------------------*/
  218.  
  219.     if (r250_index >= 147)
  220.  
  221.         j = r250_index - 147;      /* Wrap pointer around */
  222.  
  223.     else
  224.  
  225.         j = r250_index + 103;
  226.  
  227.  
  228.  
  229.     new_rand = r250_buffer[r250_index] ^= r250_buffer[j];
  230.  
  231.  
  232.  
  233.     if (r250_index >= 249)      /* Increment pointer for next time */
  234.  
  235.         r250_index = 0;
  236.  
  237.     else
  238.  
  239.         r250_index++;
  240.  
  241.  
  242.  
  243.     return new_rand;
  244.  
  245. }
  246.  
  247.  
  248.  
  249. /**** Function: r250n Description: returns a random unsigned integer k
  250.  
  251.                     uniformly distributed in the interval 0 <= k < n ****/
  252.  
  253. unsigned int r250n(unsigned n)
  254.  
  255. {
  256.  
  257. /*---------------------------------------------------------------------------*/
  258.  
  259.     register int    j;
  260.  
  261.     register unsigned int new_rand, limit;
  262.  
  263. /*---------------------------------------------------------------------------*/
  264.  
  265.     limit = (65535U/n)*n;
  266.  
  267.     do 
  268.  
  269.         {new_rand = r250();
  270.  
  271.         if (r250_index >= 147)
  272.  
  273.             j = r250_index - 147;      /* Wrap pointer around */
  274.  
  275.         else
  276.  
  277.             j = r250_index + 103;
  278.  
  279.  
  280.  
  281.         new_rand = r250_buffer[r250_index] ^= r250_buffer[j];
  282.  
  283.  
  284.  
  285.         if (r250_index >= 249)      /* Increment pointer for next time */
  286.  
  287.             r250_index = 0;
  288.  
  289.         else
  290.  
  291.             r250_index++;
  292.  
  293.         } while(new_rand >= limit);
  294.  
  295.     return new_rand%n;
  296.  
  297. }
  298.  
  299.  
  300.  
  301. /**** Function: dr250 
  302.  
  303.         Description: returns a random double z in range 0 <= z < 1.  ****/
  304.  
  305. double dr250()
  306.  
  307. {
  308.  
  309. /*---------------------------------------------------------------------------*/
  310.  
  311.     register int    j;
  312.  
  313.     register unsigned int new_rand;
  314.  
  315. /*---------------------------------------------------------------------------*/
  316.  
  317.     if (r250_index >= 147)
  318.  
  319.         j = r250_index - 147;     /* Wrap pointer around */
  320.  
  321.     else
  322.  
  323.         j = r250_index + 103;
  324.  
  325.  
  326.  
  327.     new_rand = r250_buffer[r250_index] ^= r250_buffer[j];
  328.  
  329.  
  330.  
  331.     if (r250_index >= 249)      /* Increment pointer for next time */
  332.  
  333.         r250_index = 0;
  334.  
  335.     else
  336.  
  337.         r250_index++;
  338.  
  339.     return new_rand / 65536.;   /* Return a number in [0.0 to 1.0) */
  340.  
  341. }
  342.  
  343.  
  344.  
  345.  
  346.  
  347. /***  linear congruent pseudorandom number generator for initialization ***/
  348.  
  349.  
  350.  
  351. static unsigned long seed=1;
  352.  
  353.  
  354.  
  355.     /*    Return a pseudorandom number in the interval 0 <= n < 32768.
  356.  
  357.         This produces the following sequence of pseudorandom 
  358.  
  359.         numbers:
  360.  
  361.         346, 130, 10982, 1090...  (9996 numbers skipped) ...23369,
  362.  
  363.         2020, 5703, 12762, 10828, 16252, 28648, 27041, 23444, 6604... 
  364.  
  365.     */ 
  366.  
  367. static unsigned myrand()
  368.  
  369. {
  370.  
  371.     seed = seed*0x015a4e35L + 1;
  372.  
  373.     return (seed>>16)&0x7fff;
  374.  
  375. }
  376.  
  377.  
  378.  
  379.     /*    Initialize above linear congruent pseudorandom number generator */
  380.  
  381. static void mysrand(unsigned newseed)
  382.  
  383. {    seed = newseed;
  384.  
  385. }
  386.  
  387.  
  388.  
  389.